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ABSTRACT 

Jang-Hyuk Kwon, Master of Science, 1977 

Major: Aerospace Engineering, Department of Aerophysics and Aerospace 

Engineering 

Title of Thesis: Numerical Solution of the Compressible Navier-Stokes 

Equations Using Density Gradients as Additional 
Dependent Variables 

Directed by: Dr. Joe F. Thompson 

Pages in Thesis: 75 Words in Abstract: 217 


ABSTOACT 

Numerical solution of two-dimensional, time-dependent, compress- 
ible viscous Navier-Stokes equations about arbitrary bodies is treated 
using density gradients as additional dependent variables. Thus, 

there are six dependent variables, p, u, v, E , p and p , to be com- 

s X y 

puted with the SOR iteration method. Besides the new formulation for 
pressure gradient terms, a new formulation for computing the body 
density is presented. To approximate the governing equations, an 
implicit finite difference method is employed. The coordinate system 
used here is the automatically generated body-fitted coordinate sys- 
tem that was developed at Mississippi State University. 

In computing the solution for the flow about a circular cylinder, 
a problem arose near the wall at both stagnation points. Thus, 
computations with various conditions were tried to examine the prob- 
lem. Also, computations with and without new formulations are com- 
pared. The flow variables are computed on 37 by 40 field first, then 
on an 81 by 40 field. 




V 



As a result, density profiles arc shown at different time steps 
with various conditions. Profiles for velocity and total energy, 
velocity vectors, Mach number contours and isobars are shown for 
particular cases. 

Lastly, convective terms in transformed plane due to the coor- 
dinate stretching and the transformed velocity are Introduced to 
conclude that the coordinate stretching plays an important role in 
this problem. Two ways of overcoming this problem are suggested. 
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Time step index 
Initial value 
Degrees 


Prefix 
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I. INTRODUCTION 


I An implicit finite difference method is used to solve the two- 

I 

{ • 

1 dimensional, time-dependent, compressible viscous Navier-Stokes 

[ equations about arbitrary bodies using the body-fitted coordinate 

i system. A new formulation for pressure gradient terms and body 

i density is developed and solutions are compared with those with con- 

ventional formulations. Various computations are made using differ- 
ent conditions and different ways to examine the problem. The new 
formulation is developed because the wiggles which appe^.r in the 
! results with conventional formulation can be suppressed by using 

1 the density gradient method. 

f 

i Usually, the compressible viscous flows have been solved by an 

;■ . epcpllcit scheme, and an Implicit sclieme is rarely used. Moreover, 

! 

; the compressible viscous flows about a circular cylinder have been 

rarely solved. Recent works about the flows past a circular cylinder 
; are these: [1] used polar coordinates to solve the steady incompressi- 

ble Newtonian flows with low Reynolds numbers; in [2], Impulsively 
started, time-dependent inviscid transonic flows were solved with 
different meshes, different freestream Mach numbers and Reynolds 
numbers; in [3], inviscid compressible supercritical flows were consid- 
ered using polar coordinates. The explicit scheme was used in all 
these \Torks. 


II. BODY-FITTED COORDINATE SYSTEM 



The basic idea of the body-fitted coordinate system is to generate 
the coordinate system having some coordinate line coincident with the 
boundary or body surface. A general method of generating body-fitted 
coordinate systems is to let the curvilinear coordinates to be solu- 
tions of an elliptic partial differential system in the physical planer 
with the Dirichlet boundary conditions on all boundaries. Thus, all 
the computations can be done on a rectangular fieJd with a square mesh 
and all boundary conditions can be expressed at grid points, regard- 
less of the body shape. 

Major advantages of using body-fitted coordinates are; first, the 
computer software utilized to generate the body-fitted coordinate 
system is independent of the set of partial differential equations to 
be solved on this system; second, the computer software generated to 
approximate the solution of a given set of partial differential equa- 
tions is completely independent of the physical geometry of the prob- 
lem; finally, physical Integral conservation relations need not be 
lost in the transformed plane [ 4 ], This technique is extendable to 
three dimensions and applicable to fields with time-dependent bound- 
aries. 

Consider the transformation of a two-dimensional, doubly-connected 

region, R, bounded by two closed boundaries Fj and F2 onto a rectangu- 
* 

lar region, R as shown in Figure 1 . The boundaries F^ and F2 in the 

* * 

physical plane are denoted by Fj and F2 in the transformed rectangular 
plane. 

Let us consider taking Laplace's equation with Inhomogeneous 


terms on the right hand sides as the generating elliptic system. Tluii, 


these equations are: 


'xx Sy * • 


(2.1a) 


n + n “ Q(C.n) 
XX yy ^ 


(2.1b) 


with the Dirichlet bot»ndary conditions 



m 

Cj(x,y) 

n 


- 



I- "1 


“1 



C2(x,y) 


L "2 


(x,y) e r 


1 


, (x,y) e r 


(2.1c) 


(2. Id) 


One coordinate is set to be constant on the body and the outer boundary, 
while the other is set to vary monotonically around the body. 

In the transformed plane where all computations are made, the 
transformed equations are 


ax^^ - 2Bx^^ + Yx^^ «= - j2 [x^P(C.n) + x^Q(C.n)] (2.2a) 


- 23y^^ + YX^^ “ - j2 [y^P(^,n) + y^Q(5,n)] 


(2.2b) 
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where 



with the boundary conditions 


(2.2c) 
(2. 2d) 
(2.2e) 






The functions fjCC.Hj), f 2 (^.nj), gj(C.n 2 ) and are specified 

by the known shape of the boundaries and F^, and the distribution 

of C specified. Boundary data are not needed along the re-entrant 

boundaries F^ and F, . The functions, P(C.n) and Q(C,n) are used for 
3 4 

controlling the coordinate spacing. 

The quasi-linear elliptic set of partial differential equations 
(2.2) are then solved on the rectangular transformed plane using 
finite differences and SOR iteration. 

A detailed presentation of the coordinate generation method has 


been given in [4, 5, 6, 7]. 


III. APPLICATION TO THE NAVIER-STOKES EQUATIONS 
IN PRIMITIVE VARIABLES 


The body-fitted coordinate system Is used to obtain numerical 
solutions of the twc-dimensional, time-dependent, compressible vis- 
cous Navier-Stokes equations about arbitrary bodies. The convention- 
al form of the governing equations from [8] is given first, then the 
new formulations for pressure gradients and body density are present- 
ed. In all equations, conservative forms are used. 


A. Flow Equation Formulation — Conventional Form [8] 

The primitive variable formulation of the two-dimensional, 
unsteady compressible viscous flow equations in a non-dimensional 
x,y coordinate system is given by 


. 9(pu) . 8(Pv) _ 0 
8t 3x 9y ~ 


(3.1a) 


9(pu) ^ 8(pu^) ^ 8(puv) ^ ^°xx ^ ^^yx 
3 t 3 X 3 y 3x 3y 


(3.1b) 


3(pv) . 3(pvu) , 3(pv^) ^ ^^xy ^%y 

3 t 3 X 3 y 3x 3y 


3E 
s 

3t 


3(E u) 


+ 


s 

3 X 


+ 


3(E u) 
s 

3 y 


(0 - 1)m2 


u + T v) + -;^(t 
3x XX xy 3y yx 


u + a v) 

yy J 


( ® ^ 

L (u 

+ -5- 


lPr«Ry 3 

X \ 3x / 

3y 

V 8y/. 


(3. Id) 


where 


E « p[e + (0 - 1 )m2(v 2/2)] 
s “ 


(3.2) 


b 


0, M and Pr are the ratio of heat, the freestream Mach number and 
00 

the Prandtl number, respectively. The lengths or distances are non- 
dimensional with repect to the body chord, c; time is non-dimensional 
with repect to the time required for passing one chord distance with 
the freestream velocity; density is non-dimensional with respect to 
the freestream density, p^; the velocity component u, v are non-dimen- 
sional with respect to the magnitude of the freestream velocity, V^; 
the thermodynamic verlable, the internal energy, e, is non-dimension- 
al with respect to the static enthalpy of the freestream, h^; and the 
thermodynamic variables of pressure and the total energy, E^, are non- 
dimensional with respect to the product of the freestream density and 
the static enthalpy of the freestream, P^^h^. The Prandtl number, Pr, 
is the product of viscosity and specific heat at constant pressure 
divided by thermal conductivity. Components of stress tensor are non- 
dimensional with respect to p and are given by 

00 CO 



(3.3a) 


(3.3b) 


(3.3c) 


with 


X E y* - 2y/3 and R E p V c/y 

00 00 00 

The Reynolds number, R, is based on freestream conditions and the 
body chord, c. Bulk viscosity, y' is approximately zero in the 


case of local thermodynamic equilibrium which is assumed to be the 
case in the present investigation. 

Additional relations are needed to solve the compressible flow 
equations. They are the equations of state and a relation for the 
viscosity. 


p *= (0 - l)pe 

(3.4) 

y ■ (h)^^^l^(l + hj)/(h + h^)J 

(3.5) 

In the Sutherland viscosity law, (3.5), h^ is a constant, whose value 
depends on the type of gas, and the non-dimensional enthalpy i h, 

equals the product of specific heat ratio, 0 , and 

internal energy, e. 

Boundary Conditions 


The boundary conditions on the body surface are; 

u = V = 0 

(3.6a) 

and 


T = constant 

(3.6b) 

w 


or 


q^ = 0 (adiabatic wall) 

(3.6c) 

The wall temperature, T^, is non-dimensional with 

repect to the stag- 


nation temperature of the freestream. 

Far from the body, the boundary conditions are 


P “ 1 


at infinity 


(3.7a) 




.s 


u * 1 , V “ 0 


at Infinity 


(3.7b) 


- i + (e - 


s 0 


at infinity (3.7c) 


Transformed Equations 

After transforming to the body-fitted coordinates (5,n), using 
the transformation relations of Appendix A of [4], the governing 
equations (3.1) can be written in the vector form [8]. 


3t 8n 


(3.8) 


The ordered arrays U, F, and G are: 


(3.9a) 


puu + T X - a y 

yx n XX n 

pvu + o X - T y 

yy n xy n 

Eu + Dt X - a y u + Do x 
s yx n XX n yy n 


r 

pv^ 

puv + o y_ - T X_ 

XX 5 yx C 

pvv + T y_ - cr x_ 
xy C yy ? 

Ev + Da y'.--T x^u + Dt 
^ s XX C yx C 


T y V - A e y 
xy n X n 


(3.9b) 


xy^C 


a x_ V - A e x_ 

yy C y 5 


(3.9c) 


D » (e - 1)M2 


A = 

Pr.R 


(3.10) 
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Transformed velocity components u, are defined below, and J Is the 
Jacobian of the transformation : J “ ’ 


u “ uy - vx 

n n 


V “ vx^ - uy^ 


The transformed forms of stress terms are given In equations 
(3.12), (3.13), (3.14), respectively. 


a y 
XX n 


T X 

yx C 


T y 

xy^n 


a x 

yy 5 


D -R[vc-=13“n-‘^16''c-"'l8%] 

‘’''c II r 1 

- -r R L‘=13“5 ■ ‘^ll“n - ‘=17''£ ^ ■=15%J 

px r 1 

■ ~D^ R [‘^16“5 ■ ‘^u''nj 

D * R h8“{ - ^ "12% ] 


(t X 
yx n 


a y)u+(o X -T y)v=^"7 I (c u - c ,v)u 
XX n yy n xy n d r l 9 lo £ 


+ (c V - c u)u + (c u - c v)v 

9 13 n 18 14 n 


<%0'' - "l6“>''£] 


(a y^ - T x_)u + (t y_ - o x )v “ ^ ^ (c _u - c _v)u 

xx"^C yx C xy C yy C 15 R L 13 18 £ 


+ (c,jV - Cj,U)u^ + 


+ (Cj^v - c,,u)Vj] 


e y -ex = c_e_ - c_e 
x-^n yn 25 3n 
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The subscripted coefficients c are defined in Table 3-1 with 
the definitions of a, 8, and y. Those coefficients missing from the 
table, e. g., Cj^, arise in the differencing of the transformed equa- 
tions. 


TABLE 3-1. Coefficients For Transformed Equations 


C-2 “ 

^14 “ (38 + x^x^)/3J 

C 3 » B/J 

'15 ■ Vc'” 

=4 “ y/J 

'16 ■ 

C 9 = (3a + y2)/3J 

'l7 • <*5^ ■ 

'10 ' 

'18 ■ 

Cfi “ (3 y + y|)/3J 

a = x^ + y^ 
n n 

C 12 “ (3y + x|)/3j 


Cf 3 “ (38 + y^y^)/3J 

Y - + y| 


Solution of the equations (3.8) through (3.14) uses the equation 
of state given by (3.4), the viscosity law given by (3.5) and the 
boundary conditions given by (3.6) and (3.7). There is no change in 
these equations for flow solutions in the body-fitted coordinate 
system. 

B. New Formulation of Flow Equations 

Generally, the pressure gradient terms in transformed momentum 
equations are computed by using the equation of state for pressure 
values and central difference approximations for derivatives. 
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New formulations for computing the pressure gradient terms and the 
body density are discussed below. 

Formulation of Equations for Pressure Gradients 

The equations of state (3,4) is used and differentiated as a 
whole with respect to x and y to get the pressure gradient terms. 
Thus, the pressure gradient terms of x and y momentum equations in 
(3.1) and (3.3) can be replaced by 


i£ 

3x 


(0 - 1) [P 


3x 




(3.15a) 


l£ 

3y 


(0 - 


1 \ r 3® 


3y' 


(3.15b) 


The transformed derivatives of internal energy are expressed as 
below. The Internal energy at each point can be calculated from the 
energy equation. 


3e 1 

^ “ 1 f 


(3.16a) 


i®. 

3y 





X e_l 

n 


(3.16b) 


To get p and p in (3.15), two equations are obtained from the 
X y 

continuity equation by differentiation with respect to x and y, re- 


spectively: 
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at 


a^(pu) a^(pv) 

3x2 ' 


(3.17a) 




(3.17b) 


After transformation, the following expressions for p^ 


and p 

y 


are 


obtained. 


3p 

= - RUDXX -RVDXY 

a t 


“ - RUDYX -RVDYY 


where 



(3.18a) 


(3.18b) 


(3.18c) 


(3.18d) 


RUDYX = 


’'A T Kn 

t) 


RVDYY 


-{t)A|t<»>, -?<•■’. 


- (i) 


(3.18c) 


(3.18f) 


Thus, p and p can be obtained from the finite difference ver- 
X y 

sion of the following equations. 


(P - (P )"■ 

X 1, i X i, 


- RUDXXI" - RVDXyI" (3.19a) 

i » j i > J 




= - RUDYXl" - RVDYyI" (3.19b) 

i > j • i > J 


After dropping the superscript n and denoting the superscript (n-1) 
by o. 


(p ). . 
X i,j 


(Px>° j - At [^RUDXXl^ j + RVDXyI^ ^ j (3.20a) 


Vl.j ■ 


(3.20b) 


p and p are treated as new dependent variables which change with 
X y 
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time. They are solved with the same procedure as the other dependcuit 

variables, p, u, v and E , using the SOR iteration method. 

s 

Thus, the pressure gradients, p and p , in (3.15) can be obtained 

X y 

directly using the equations (3.16) and (3.20) for tlie computations 

of e , e , p and p , rather than by differencing the pressure. Then, 
X y X y 

the momentum equations in (3,1) are changed as follows : 


3<pu) . 

3(pu^) 

a(puv) _ 

aa' ax 

1 ap , XX , vx 

a t 

ax ^ 

a y 

■ D ax ax Ty~ 

3(pv) 1 

3(pvu) 

a(pv^) ^ 

, . 3t ao' 

_ 1 

a t 

ax 

a V 

D ay a X ay 


(3.21a) 


(3.21b) 


where 


o' 

XX 


a' 

yy 


A + Iv \ 4 . Ji /oIh \ 
R \3x 9y y R rax / 

+ 9v\ £ / ^\ 

R Vax ay / R Vdy ) 


(3.22a) 


(3.22b) 


The rest are all same as (3.1), (3.2) and (3.3). In transformed equa- 
tions, (3.15) is used for pressure gradients, and o' and o' are used 

XX yy 

instead of a and a 

XX yy 

Other relations and boundary conditions are the same. When iter- 
ating for u and v, the terms p and p are computed using the most 

X y 

recent iterate values of p , p and e , e . 

X y X y 


Formulation of Equations for Body Density 


1 


To compute the body density, either Llie t.ontimiily eqii.it ion , or 
momentum equations evaluated on the body, or some type of the extra- 
polation method has commonly been used. Here, another way of comput- 


ing the body density is presented by using the values of p and p 

X y 

computed above from (3.20). 

In the transformed plane, p and p are expressed as 

X y 


Px " jlVc’Vnl 


(3.23a) 


p - T [-X p_ + x_p ] 

y J ^ n C C n' 


(3.23b) 


From these equations, (3.23a) and (3.23b), multiplied by x^ and y^ 
respectively, the following relation can be derived. 


X p + y p 

n X n y 


(3.24) 


Thus, p is expressed by a function of p and p . Now, p can be 
n ' X y n 

expressed as below using the second order central difference approx- 
imation. 


n'i,2 


2 t^i,3 " ^i,l^ 


(3.25) 


This gives the body density in the form of 


*i,l “ ^i,3 " ^ %'i,2 


(3.26a) 


or. 


(3.2(ib) 




1.1 


p -2(xp + yp||,o 

1,3 ri X ■'ti y' ' i ,2 


Equation (3.26b) is used for the body density with 2 ‘^^Iculated 

by the equation (3.2A) with the most recent iterate values of (p ^ 

and (p ) Again, the 50R technique is used to compute the body 

y i, 2 

density. 

In the new formulation of flow equations, the pressure gradient 
terras are not represented by pressure differences, but are evaluated 
directly by using the values of energy gradients and density gradi- 
ents. 




•1 
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IV. NUMERICAL rOR.'n;iJVTTON 

All space derivatives in the field are approximated by the 
second-order, central difference expressions. Thus, the first and 
second derivatives are approximated as below setting and An to be 
unity. 


(f ) * — (f - f ) 


^^n^l,j “ 2 ^^l,j+l " ^i.j-1^ 


(f ) * f - 2f + f 

^ CC^l.j i+lj i.j 1-1, J 


* ^i,j+l " ^^i,j ^l,J-i 


(A.l) 

(4.2) 

(4.3) 

(4.4) 


* 4 ^^1+1, j+1 " ^i+l,j-l " ‘'i-l.j-P 


At the body surface, where central difference approximations can not 
be used for the derivatives with respect to n, the second- order, 
one-sided difference expressions are used. 




(4.6) 


In the ^[-directions, the central difference expressions are also used. 

All the time derivatives are approximated by the first-order, 
backward difference expressions to get the implicit formulation. 


ii. 
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(f 


n _ n-1 

JLJ 

A t 


(4.7) 


The explicit method was not used because with explicit method, 
the elliptic character of spatial variation is totally neglected. 
Moreover, explicit methods have severe stability limitations [5]. 

Thus, the implicit formulation is used to solve the Navler- 
Stokes equations, and the set of resulting simultaneous difference 
equations with the boundary conditions are solved at each time step by 
the point SOR iteration. 

For convenience sake, the finite difference equations for pres- 
sure gradient terms and body density are explained briefly in Chapter 
III. All other equations in difference expressions are as in [8]. 


V. RESULTS AND DISCUSSION 


f 


t 


I 


[ 


f 

i 


ly 


The case considered here Is the flow about a circular cylinder 
at a freesream Reynolds number of 1000. The Mach number and the 
static temperature of the freestream were 0.8 and 273“K, respectively. 

The Prandtl number was set at 0.71 [9], and the viscosity law refer- 
ence temperature, hj^ in (3.5), was 0.40293 [10], For an initial 
solution, the incompressible potential flow solution was used. The 
field acceleration parameters were varied, but the convergence toler- 
ance for all variables were set at 0.00001. The time steps used here 
were all 0.01, except in the runs of 3. 

Runs were made on two different size fields; 37 by 40 field and 
81 by 40 field, so the number of fi-llnes was increased later as sho\vm 
in Figure 2 and 3. The field radius was 6.0 in either field. Log 
attraction was used to locate 25 lines inside the R=20 boundary layer 
in Figure 2 and 17 lines inside the R=106 boundary layer :’n Figure 3 
by using [11]. 

As soon as the first run with the new formulation for pressure 
gradient terms and the continuity equation for the body density was 
made, dips and spikes showed up near the wall and developed as can be 
seen in Figure 4 and 5. 

All the computations except 6, were made with the new formula- 
tion for pressure gradient terms and all the computations, except 1 
and 6, were made by using the new formulation for body density. The 
adiabatic wall was assumed in all cases except run 4. 
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A. Computations Made with Various Conditions 
On 37 by 40 Field 

1. Computation with the Continuity Equation for Body Density 

At first, the new formulation was substituted for the pressure 

gradient terms in the conventional formulation. The continuity equa- 
tion evaluated at the wall was used for the body density as in [8]. 

The acceleration parameters for density, velocity, total energy, 
and were 0.95, 0.50, 0.80, 0.40 and 0.50, respectively. The accel- 
eration parameters were set at 0.30 for body density and 1.0 for free- 
stream variables. 

As soon as two time steps ( t = 0.02 ) passed, dips and spikes 
appeared in density profile and developed as the time proceeded. The 
density profiles are sho\im in Figure 4 and 5. Near the front stagna- 
tion point (C “ 28 line), the problem seems to be more serious than 
near the rear stagnation point (C = 10 line). 

2. Computation with the New Formulation for Body Density 

The new formulation for body density was tried from this compu- 
tation to the end, except the computation 6. Other conditions were 
the same as those of run 1. The density profiles are shown in Figure 
6 and 7. Runs were made up to 1.0 in non-dimensional time. The pro- 
files look better than those of run 1. The dips near the rear stag- 
nation point became smaller as the time passed, but not small enough. 
Near the front stagnation point, the spikes developed and the differ- 
ences of the values between the wall and the next point remained 


almost the same. 
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3. Computations) with Different Time Steps 

Time steps of 0.02 and 0.005 were tried using same conditions of 
run 2. As shown in Figure 8, 9 and 10, the results were almost same 
as that of run 2. The acceleration parameter of 1.0 was used for body 
density. 

4. Computation with a Constant Wall Temperature 

Instead of the adiabatic wall, the constant temperature wall was 
assumed. The time step was 0.01 and still the new formulations for 
pressure gradient terms and body density were used. Other conditions 
were the same as those of run 3. Wall temperatures were fixed at 1.0. 
As can be seen in Figure 11, the results became even worse. The values 
of body density went down at both stagnation points. 

On 81 by 40 Field 

5. Computation with New Formulations 

The coordinate system was changed so as to. get smaller grid sizes 

in C direction. The time step was 0.01 and acceleration parameters for 

density, veloslty, total energy, p , p and body density were 0.95, 

X y 

0.70, 0.80, 0.70, 0.60 and 1.0, respectively. The density profiles are 
shown in Figure 12 and 13. As can be seen from the figures, the dips 
were even worse than those of run 2 which was made on the 37 by 40 
field with the same conditions. Still there are no oscillation except 
the dips and spikes. As the time passes, the dips remain almost same 
near the front stagnation point ( C *= 61 line) and slightly Increasing 
near the rear stagnation point ( C = 21 line). 

• In Figure 14 and 15, velocity profiles and total energy profiles 
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are shown. From the figures, it can be observed tliat the velocity 
profiles near the stagnation points are smooth and have no dips while 
the other have dips. Also, in Figure 16, the velocity vectors are 
shown near the front stagnation point. Applying the continuity equa- 
tion in difference form to a point (61,2), lying next to the front 
stagnation point to n-direction, then it is easy to see that the out- 
flow to n-direction is much more than the inflow, so that makes the 
density at this point go down. 

6. Computation with Conventional Formulations [8] 

For making a comparison with new formulations, the formulation of 
[8] was run at the condition of run 5. As can be seen In Figure 17 
and 18, the density profiles look better near the wall than those of 
run 5 and other portions of density profiles are same. The oscilla- 
tions near the wall occur in this computation, but not in the computa- 
tion with the new formulations. When time is 1.0, the oscillations 
still exist near the wall though they are small. They are shown up to 
more than five cell width, sometimes more than ten cell width. But 
the sizes of the oscillations can be considered as a constant. 

The computer time required was slightly more than a half of run 
5. The acceleration parameters for density, velocity, total energy 
and body density used were 0.95, 0.70, 0.70 and 0.50, respectively. 

The velocity vector plots with new formulations on 37 by 40 field 
and 81 by 40 field are shown in Figure 19 and 20, respectively. In 
Figure 21, the velocity vector plots with conventional formulations are 
shown. In the rear region, the formulations of circulating eddies can 
be seen. In the plot with new formulations on 81 by 40 field, waves 
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appeared. They appeared to be related to shock waves and velocity 
vectors oscillate near the waves. In Figure 22 and 23, Mach number 
contours and Isobars are shown for the waves. From a study of the 
Mach number contours, it looks that shock waves will later be formed. 
The Mach number downstream of the waves are too high for these waves 
to be true shocks. Isobars are so close to each other near the waves. 
In plots with new formulations on 37 by AO field and with conventional 
formulations on 81 by AO field, the waves are smeared out. 

7. Computation with Reynolds Number of 50 

The density profiles are shown in Figure 2A. New formulations 
for pressure gradient terms and body density were used with same con- 
ditions as run 5. This computation was made up to time of 0.1 to see 
whether the viscosity can cure the problem or not. As can be seen in 
the figure, dips and spikes are smoothed out. 

B. Convection Terms Check 

Besides the runs with various conditions, the convection terms 
were examined. At first, the expanded forms of convection terms in 
X and y momentum equations were used to check which terms are changed 
by the expanded forms. The expanded form used here is, for example, 
as below. 

3 (uv) 3v , 3u 

3 X 3x 3x 

All the terms that were calculated with and without using the expand- 
ed forms were computed. As a result, the changes in convection terms 
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near the stagnation points were pretty big. Thus, two changes below 
were tried to correct the problem. 


1. Averages In the Continuity Equation 

The continuity equation in Integral form Is 


// p^J d^dn “ - / pV-dA 


(5.1) 


where y is the transformed 

A ^ A 

velocity expressed by V = ui + v^ . 

The integral equation can be 
approximated as below: 
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In this way, p^ was calculated to get p^, but the result was even 
worse and densities went down more. 

2. POA Forms 

The product of average forms were tried in the continuity equa- 
tion and X and y momentum equations and used only at the front and i 
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rear stagnation point because the POA form can approximate better 
in the flow reversal region. But there were only very small changes 
in flow variables or no change at all. 

C. Cu , Cv and Ce 

n n n 

From the results of various computations and the convection terms 
check, the following conclusions can be drawn : 

1. The problem of dips and spikes is not an instability problems 
because at later time the magnitudes of dips remain almost same. In- 
stability usually means oscillations and a rapid growth of errors. 

2. There are truncation errors in computations near the stagna- 
tion points, but those are not due to the difference expressions of 
derivatives because the errors remained nearly constant when the cell 
size in ^ direction was reduced. 

From [8], Cu , Cv and Ce are introduced as below by taking domi- 
nant terms in x and y momentum equations and energy equation near the 
stagnation point : 



Where c,,, c,. and c, are shown in Table 3-1. These terms act like 
11 12 4 


convective terms in the transformed plane. The first terras are due 


to the coordinate streching In n direction .ind the second terms are 




due to the transformed velocity. Because dips are in density profiles 

and energy profiles, Ce is calculated as in Table 5 - 1 and pe 

n 

profiles are shown in Figure 25. The calculations were made with new 
formulations at t ■ 1.0. 


TABLE 5-1. Cep 


I 

J 

Ce 

n 


A 

- pv 
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1 

- 2.142 

X 

lolj 

0 


-4 


2 

3.341 

X 

10.3 

5.060 

X 

10_3 


3 

1.824 

X 

10 , 

1.970 

X 

10 ^ 


4 

4.014 

X 

loi 

4.158 

X 

10-3 


5 

6.652 

X 


6.803 

X 

■o: 


6 

9.542 

X 

10 

9.702 

X 

10 

67 

1 

- 2.135 

X 

lolj 

0 


—4 


2 

2.110 

X 

10 - 

4.081 

X 

10 - 


3 

1.372 

X 

-3 

10_3 

1.548 

X 

■«: 


4 

3.053 

X 

10 :: 

3.209 

X 

10 ^ 


5 

5.051 

X 

loi 

5. 198 

X 

10-3 


6 

7.243 

X 

10-3 

7.388 

X 

10-3 

71 

1 

- 2.024 

X 


0 


-4 


2 

- 6.198 

X 

^°-4 

1.874 

X 

10_4 


3 

4. 160 

X 

10. 3 

6.574 

X 

10_3 


4 

1.037 

X 

10 :: 

1.223 

X 

10 :: 


5 

1.557 

X 

loi 

1.782 

X 

10-3 


6 

2.431 

X 

10-3 

2.560 

X 

10-3 


Near or on the body surface, the first term is dominant while in 
the other region, the second term is dominant. In [8], the region 
where the first term is dominant is wider than this case because in 
his case the Reynolds number is much higher. In Figure 23, the pe 
profiles are not realistic, too. 


From these calculations, it can be said that the coordinate 
stretching makes the problem, or, at least, plays an important role 
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in the problem. 

In regions of large gradients, a fine mesh is usually used. In 
[12, 13], halving the coarse mesh cell size was used to get a medium 
mesh and a fine mesh, and the governing equations were solved sepa- 
rately. In [14], a fine, exponentially stretched mesh spacing was 
employed and explicit artificial viscosity terras were added to avoid 
the instability. In [15, 16], the same stretched coordinate spacing 

! 

as [14] was used but a fine constant mesh spacing was used in the jj 

i 

sublayer region. I 
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VI. CONCI.USION 

An implicit finite difference method was applied to solve the 
compressible viscous flows about a circular cylinder using new formu- 
lations for pressure gradient terras and the bedy density, but a prob- 
lem arose near the stagnation points in both results with and without 
the new formulations. In the result with conventional formulation 
[8], wiggles appear near the wall, but not in the results with new 
foirmulatlons, except the dips and spikes. It was concluded that the 
coordinate stretching in the n-dlrectlon plays an Important role on 
the truncation errors which causes the problem. 

This problem can be corrected by; 1) employing explicit artifi- 
cial viscosity terms, or 2) employing a few n lines which have a 
constant or very slowly changing spacing for a first few rj lines. 
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Figure 2. Coordinate System - 37 by 40 Field 



Figure 3. Coordinate System - 81 by AO Field 
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Figure 4a. Density Profile with New Formulation for Pressure 
Gradients at t = 0.1 (37 by 40 Field) 
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(b) Along C “ 10 line 


Figure 4b. Density Profile with New Formulation for Pressure 
Gradients at t = 0.1 (37 by 40 Field) 
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(a) Along 5 = 28 line 


Figure 5a 


Density Profile with New Formulation for Pressure 
Gradients at t = 0.3 (37 by 40 Field) 
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(b) Along 5 = 10 line 

Figure 5b. Density Profile with New Formulation for Pressure 
Gradients at t - 0.3 (37 by' 40 Field). 
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(a) Along 5 = 28 line 


Figure 6a. Density Profile with New Formulation for Body 
Density at t =* 0.3 (37 by 40 Field) 
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(b) Along 5 = 10 line 


Figure 6b. Density Profile with New Formulation for Body 
Density at t = 0.3 (37 by -40 Field) 
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Figure 7a. Density Profile with New Formulation for Body 
Density at t = 1.0 (37 by 40 Field) 
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(a) Along 5 = 28 line 


Figure 8a. Density Profile with At = 0.02 at t * 0.3 (37 by 
40 Field) 
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; (b) Along 5 “ 10 line 
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• Figure 8b. Density Profile with At = 0.02 at t = 0.3 (37 by 
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(b) Along C = 10 line 


Figure 9b. Density Profile with At = 0.005 at t = 0.3 (37 by 
40 Field) 





(a) Along C = 28 line 


Density Profile with t = 0.005 at t 
by 40 Field) 
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(a) Along i - 28 line 

Hgure 11a. Density Profile with A Constant Temperature Wall 

at t = 0.3 (37 by 40 Field) - 
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(b) Along ^ 10 line 

Figure lib. Density Profile with A Constant Temperature Wall 
at t = 0.3 (37 by 40 Field) 
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Figure 12a. Density Profile with New Formulations at t - 0.3 
(81 by 40 Field) 
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(a) Along C “ 61 line 

Figure 13a. Density Profile with New Formulations at t = 1.0 
(81 by 40 Field) 
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(b) Along C = 21 line 

Figure 13b. Density Profile with New Formulations at t = 1.0 
(81 by 40 Field) 
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(a) Along C = 61 line 


Figure 14a. Velocity Profile with New Formulations at t *» 1,0 
(81 by 40 Field) 
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Figure 15b. Total Energy Profile with New Formulations at 
t = 1.0 (81 by 40 Field) 














(a) Along 5 ** 61 line 


Figure 17a. Density Profile with Conventional Formulations at 
t = 0.3 (81 by AO Field) 
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(b) Along E = 21 line 


Figure 17b. Density Profile with Conventional Formulations at 
t = 0.3 (81 by 40 Field) 
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(a) Along C = 61 line 

Figure 18a. Density Profile with Conventional Formulations at 
t = 1.0 (81 by 40 Field) 
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(b) Along 5 = 21 line 


Figure 18b. Density Profile with Conventional Formulations at 
t = 1.0 (81 by 40 Field) 
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(a) t = 0.3 


Figure 19a, Velocity Vectors with New Formulations on 37 by 
40 Field 
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(a) t = 0.3 


Figure 20a. Velocity Vectors with New Formulations on 81 by 
40 Field 
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(b) t = 1.0 

Figure 20b. Velocity Vectors with New Formulations on 81 by 
40 Field 
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(b) t = 1.0 

Figure 21b. Velocity Vectors witli Conventional Formulations 
on 81 by 40 Field 
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Figure 22b. Mach Number Contours with New Formulations on 81 
by 40 Field 
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Figure 23. Isobars with New Formulations on 81 by 40 Field at t 
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(a) Along C == 61 line 

Figure 24a. Density Profile with R = 50 at t = 0.1 (81 by 40 Field) 
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